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We derive semiclassical transport equations for a trapped atomic Fermi gas in the BCS phase at 
temperatures between zero and the superfluid transition temperature. These equations interpolate 
between the two well-known limiting cases of superfluid hydrodynamics at zero temperature and the 
Vlasov equation at the critical one. The linearized version of these equations, valid for small devia- 
tions from equilibrium, is worked out and applied to two simple examples where analytical solutions 
can be found: a sound wave in a uniform medium and the quadrupole excitation in a spherical har- 
monic trap. In spite of some simplifying approximations, the main qualitative results of quantum 
mechanical calculations are reproduced, which are the different frequencies of the quadrupole mode 
at zero and the critical temperature and strong Landau damping at intermediate temperatures. 
In addition we suggest a numerical method for solving the semiclassical equations without further 
approximations. 
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I. INTRODUCTION 

Due to improved cooling techniques, current experi- 
ments with trapped fermionic atoms like 6 Li or 40 K reach 
very low temperatures of the order of T 0.03Tp p|, 
where Tp = ep/ks denotes the degeneracy tempera- 
ture. The main motivation for these experiments is to 
study the so-called BEC-BCS crossover by tuning the 
magnetic field around a Feshbach resonance, thus chang- 
ing the atom-atom scattering length a from the repul- 
sive side (a > 0) through the unitary limit (a — > oo) to 
the attractive side (a < 0). On the BEC side, where 
the system forms a Bose-Einstein condensate (BEC) of 
tightly bound molecules, as well as on the BCS side of 
the crossover, where the atoms form Cooper pairs which 
are very large compared with the mean distance between 
atoms, one expects that the system is superfluid, pro- 
vided the temperature lies below a certain critical tem- 
perature T c . Until now the experiments have concen- 
trated on the crossover region, but the low temperatures 
which are currently reached suggest that future experi- 
ments will also be able to study the superfluid BCS phase 
(a < and kp\a\ <C 1), although its critical temperature 
will be extremely low. 

In order to find signals for superfluidity, some re- 
cent experiments with ultracold trapped Fermionic atoms 
looked at dynamical observables like the expansion of the 
atom cloud after the trap has been switched off 0, or 
collective oscillations of the cloud Q, 0] • The theoretical 
interpretation of such experiments is usually based on a 
theory called "superfluid hydrodynamics" [J, [a, [(j , which 
is valid for a superfluid at zero temperature if the trap 
potential is sufficiently smooth to justify local-density ap- 
proximation. Apart from the latter condition, which is 
not necessarily fulfilled in the experiments 0, it is also 
clear that experiments are not done at zero temperature. 

Very recently, Landau's two-fluid hydrodynamics has 
been used to describe collective modes in trapped super- 
fluid gases at non-zero temperature @ . In this approach, 
in the temperature range < T < T c , a certain fraction 



of the atoms is not superfluid but forms a normal-fluid 
component with density p n , while the remaining atoms 
with density p s = p — p n still behave like a superfluid. 
It is also assumed that the atoms undergo enough colli- 
sions to be always in local equilibrium. This condition, 
however, cannot be taken for granted. In Ref. |9| it was 
found that even in the unitary limit, where the scattering 
length a diverges, the collision rate might be too low to 
ensure hydrodynamic behavior of the normal phase. This 
is certainly true in the BCS phase, where kp\a\ and the 
temperature are so small that one can safely assume that 
the system is in the so-called collisionless regime. Colli- 
sionless means in this context that the collision rate 1/r 
is much smaller than the trap frequency f2, i.e., an atom 
performs several oscillations in the trap before colliding 
with another atom. Since the frequencies of the collective 
oscillations are of the order of the trap frequency f2, this 
implies that it is impossible to reach local equilibrium 
during the oscillation. 

Nevertheless the idea of a two-fluid model is use- 
fid in the collisionless regime, too. It has been devel- 
oped for this case in the theory of superconductivity 
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Il3| . Similar approaches to describe liquid 
3 He should be mentioned as well, although they are more 
complicated because of the spin structure of the order pa- 
rameter 113 ■ Recently the two-fluid model has also 
been applied to the case of trapped fermionic atoms in 
the BCS phase 0, 0|. Because of the possibility of 
Fermi-surface deformations, the normal component of a 
collisionless gas docs not behave hydrodynamically, but 
more like an elastic body. A semiclassical method for 
treating the Fermi surface deformation in a normal Fermi 
gas is given by the Vlasov equation. The latter was used 
with great success in nuclear physics, e.g., in order to 
describe giant resonances in atomic nuclei [Isj . and re- 
cently it was also applied to trapped atomic Fermi gases 
in order to predict the frequencies of collective modes in 
the collisionless regime £J . Contrary to hydrodynamical 
equations, where all quantities are local (i.e., functions 
of the spatial coordinate r only), the Vlasov equation re- 
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quires a phase-space description (i.e., the quantities are 
functions of r and p). The aim of the present article is to 
derive a hydrodynamical equation for the superfiuid com- 
ponent coupled to a Vlasov equation for the normal com- 
ponent, interpolating between superfiuid hydrodynamics 
at T = and the usual Vlasov equation at T = T c . In 
principle, as it was done in the theory of liquid 3 He |l4j . 
one could also think of including a collision term into 
this equation, in order to treat systems which are nei- 
ther collisionless nor hydrodynamical, but somewhere in 
between. However, in the present article we will restrict 
ourselves to the collisionless case. 

Like the semiclassical description of the ground state 
(Thomas- Fermi approximation), the semiclassical de- 
scription of the dynamics of the system can be expected 
to become more and more accurate if the number of 
atoms in the trap increases. This was the main moti- 
vation for us to develop the semiclassical approach pre- 
sented here. A fully quantum-mechanical description 
of the collective modes of a trapped Fermi gas can be 
obtained, e.g., by the quasiparticle random-phase ap- 
proximation (QRPA), corresponding to the lineariza- 
tion of the time-dependent Bogoliubov-de Gennes equa- 
tions around equilibrium. The latter are also known 
as time-dependent Hartree-Fock-Bogoliubov (TDHFB) 
equations, especially in nuclear physics. QRPA calcula- 
tions become tremendously difficult and time-consuming 
if the number of particles increases. At present, they are 
restricted to systems of ~ 10 4 atoms 0, 0, EJ, while 
the numbers of atoms in the experiments are at least 
ten times larger. In addition, all present QRPA calcu- 
lations are done for the case of spherically symmetric 
traps, while the traps used in the experiments are gener- 
ally not spherical. The numerical solution of the QRPA 
equations without spherical symmetry seems to be almost 
unfeasible, unless one reduces drastically the number of 
particles. Therefore semiclassical approaches are at the 
moment the only way to perform calculations for large 
numbers of atoms in realistic trap geometries. 

Our article is organized as follows. In Sec. [H] we will 
present the formalism. Having derived a quasiparticle 
transport equation in Sec. Ill Al an important point will 
be to work out the linearized version of this equation in 
order to apply it to oscillations around the equilibrium 
state. This is done in Sec. IIIBI In Sec. Ill CI we will 
show explicitly that our equations indeed reproduce su- 
perfiuid hydrodynamics and the Vlasov equation in the 
limits of zero and critical temperature, respectively. The 
next part, Sec. IIIII is devoted to two simple examples 
for which our equations can be solved more or less an- 
alytically. The first example, discussed in Sec. IIII Al is 
a sound wave in a uniform gas. The second one, de- 
scribed in Sec. IIII Bl concerns a quadrupole oscillation of 
a harmonically trapped gas with some additional simpli- 
fications. Finally, in Sec. II VI we will summarize and draw 
our conclusions. 



II. FORMALISM 
A. Derivation of a quasiparticle transport equation 

In this subsection we will derive a quasiparticle trans- 
port equation for a superfiuid gas of trapped fermionic 
atoms in the BCS phase. Throughout this article we 
will assume that the two spin states f and J, are equally 
populated, which allows us to remove the spin degree of 
freedom from the beginning. However, the generalization 
to include the spin, which in fact would be necessary, e.g., 
in order to describe spin waves or systems with unequal 
populations, is straight-forward. In order to be in the 
BCS phase, the atoms must have an attractive interac- 
tion, i.e., a negative scattering length a < 0, which on 
the other hand must be weak enough for the BCS ap- 
proximation to be valid. 

Let us start by writing down the TDHFB equations 
[l8j . To that end we define the non-local normal and 
anomalous density matrices, 

g(r,r') = (^(r')^W) = tyj (r>i W) , (1) 
k(t, r') = (V> T (r>| (r)) = - (r')Vt «) , (2) 

where ip is the field operator. The single-particle hamil- 
tonian (minus the chemical potential /i) reads 

h 2 V 2 

h= = + V ext (r) +gp(r)-n, (3) 

2m 

where m is the atomic mass, V ex t(r) is the potential of 
the trap, gp(r) is the mean-field potential. The coupling 
constant g is related to the atom-atom scattering length 
a by g — Awh 2 a/m and the density per spin state p{r) is 
just equal to the local part of the density matrix, 

p(r) = g(r, r) . (4) 

According to the usual regularization prescription |2l| . 
the pairing gap is related to the anomalous density by 

A(r) = —q lim — s k(t H — , r V (5) 

v ' ' v s^q ds V 2' 2 J y ' 

Combining all quantities in the 2x2 matrices 

*=U. A 5 )- *-(-'• I-*)' (6 > 

where g and h denote the time-reversed operators to g 
and h, respectively, the TDHFB equation can be written 
in the compact form 

iKR = [H,K] . (7) 

In analogy to the derivation of the Vlasov equation in 
the normal phase from the Hartree-Fock equation , it 
is useful to introduce the Wigner transform of the density 
matrix, 

g(r, p) = J dhe-^ s ' h g(r + |, r - |) . (8) 
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It is appealing, although strictly speaking not correct, to 
interprete the function g(r, p) as a distribution function 
of particles in phase space. In a completely analogous 
way we define the Wigner transform of the anomalous 
density matrix, k(t, p), and the Wigner transform of the 
hamiltonian, h(r, p), which is equal to the classical hamil- 
tonian 



h(r,p) 



P 

2m 



+ V ext {v) + gp{v) - p. 



(9) 



(For the sake of readability we are using the same sym- 
bol for the operators and their Wigner transforms, but 
whenever there is a risk of confusion we will write down 
the arguments.) Eqs. (0J and JSJ can be written in terms 
of the Wigner transformed quantities as follows: 



d 3 p 



p{v) =] (2^ (r ' P) ' 

d 3 p ( ( _ _^ A(r) 



f d P ( , \ * 



jm 



(10) 
(11) 



We also need the Wigner transforms of the time-reversed 



operators g and h, and the Wigner transforms of the 
adjoint operators and A T . To that end we recall the 
general relations 

A(r,p) = A(r,-p), [yl+](r,p) = yl*(r ) p) ) (12) 

which are valid for an arbitrary operator A. The use- 
fulness of the Wigner transform lies in the fact that, to 
first order in an expansion into powers of h, the Wigner 
transform of the product of two operators A and B can 
be obtained according to 

[AB] (r, p) a A(r, p)B(r, p) + j{A(r, p),B(r, p)} , 



(13) 



{A B} = £ (**°° ,14-, 



where {•, •} denotes te Poisson bracket 

/dAdB 6 
\ dri dpi dpi dti ) 



Applying this product rule to the Wigner transform of 
the TDHFB equation 10, one obtains four coupled equa- 
tions: 



ihg = ih{h, g} + 2ilm(A* k) - ihRe{A* , k} , (15a) 

ihk = (h + h)n + — h,K} + A(g + g - 1) - — {A, g - g} , (15b) 

ihk* = -{h + h)n* + y {h - h, n*} -A*(g+g-l)- y {A*, g-g}, (15c) 

ihg-= -ih{h, g} +2iIm(A*re) + ihRe{A* , k} . (15d) 



In order to proceed further, it is useful to separate 
the collective superfluid motion from the dynamics due 
to quasiparticle excitations. This can be achieved by a 
gauge transformation, i/j(r) = ip(r) exp[z<^(r)] 0, ll^| . 
According to their definitions, the normal and anoma- 
lous density matrices behave very differently under this 
transformation. The corresponding Wigner transforms 
are given by 

e(r,p) = e[r,p-»V^(r)] ) (16) 
£(r,p) = K (r,p)e 2 ^ r >. (17) 

If the hamiltonian and the gap are changed according to 

~ [p_W^_^ 

2m 

A(r) = A(r)e 2i *< r >, (19) 

the equation of motion of the gauge transformed quan- 
tities looks exactly like Eq. J7J. The superfluid velocity 
is proportional to the gradient of the phase of the gap. 
Hence, if we choose the gauge transformation such that 



the transformed gap A is real, we have completely sepa- 
rated the collective motion of the superfluid component 
from the motion due to quasiparticle excitations. A for- 
mal argument for the necessity of this choice of the gauge 
is given in Ref. [T^j. 

From now on we will suppose that A is real. Splitting 
g and h into time-even and time-odd parts, 

Qev = ^{g + g), Qod = 7;(p- g) , (20) 
~K l = hh-h)=-- P -V4,, (22) 

2 v ' m 
and ii into real and imaginary parts, 

—— .Re tv , ^ivn — — Im Kj , (23) 

one can rewrite the gauge transformed version of the sys- 
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tern of equations (|15|) as follows: 

hg ev = h{h ev , Q od \ + h{h od , g ev } + 2Ak lm , (24a) 
hg od = Hh 

ev i Qev 

} + h{h od ,g od }-h{A,k re }, (24b) 

hk re = 2h ev k lm + h{h od , k re } - h{A, g od } , (24c) 

hk„ n = -2h ev k re + A(l - 2g ev ) + h{h od , k m } . (24d) 

For a semiclassical h expansion it seems disturbing that 
these equations mix different orders in h. However, it 
is possible to decouple the equations of motion for the 
leading-order quantities from those of the higher-order 
ones. In order to show this, we expand g and k into 
powers of h. Since the Eqs. (|24H themselves are only 
valid up to order h, it does not make sense to go beyond 
the first order in this series. From Eqs. (|24a|l and (|24cl) 
it is evident that ki m must be suppressed by one power 
of h with respect to the other quantities. We therefore 
write 



- ~(o) h ~(i) 

Qev,od — Q e v,od "+~ a 6ev,od ' ' ' ' > 
kre = K>rp ~t~ hk\ ) ~l~ • • • 5 



f^im. — hfZ 



a« 



(25) 
(26) 
(27) 



Inserting these expansions into Eqs. I|24a|l - i|24d|) and 
retaining only the leading order in each equation [order 
h in the case of Eqs. (|24a[) - (|24c|) . order 1 in the case of 
Eq. I|24d(l ]. one obtains 



§£> = {h ev , g^} + {h od , qM } + 2Ak^ , (28a) 

= {hv, $>} + {hod, g o h - {A, 4?} , (28b) 
k r °J = 2h ev k^ n + {h od , k r °J} - {A, g^} , (28c) 
2^ (nl /cW=A(l-2eW). (28d) 

Only one of the higher-order quantities, namely k^, ap- 
pears in these equations, but it can be expressed in terms 
of the leading-order quantities, e.g., with the help of Eq. 

«im = 7^ (e e °J - {he V , g { o h - {Kd, ~g { e °J}) . (29) 

By taking a linear combination of Eqs. (|28a|l and (|28c|) 
one can eliminate kQ. The resulting equation reads 

hevQ^ ~ ^ = E ev {E ev ,Q { °}} + h e v{h od , g e ° v ] } 

- A{h d, k r °J} , (30) 
where we have introduced the abbreviation 



A 2 



(31) 



Eqs. H28b|) . H28d.pt . and pU)l form a system of three cou- 
pled equations for the three leading-order quantities g 
g ( ° d \ and k r e '. 



From now on we will suppress the index "(0)" and 
simply write g ev , g od , and k re instead of g e °v , and 

Kre^- The next step is to exploit Eq. i|28d(l in order to 
reduce the number of unknown functions. To that end 
we introduce a new phase-space function v ev (r,p), the 
so-called "quasiparticle distribution function" , which is 
defined in such a way that the two members of Eq. (|28d|) 
are equal to h ev A(l — 2v ev )/E ev . In other words, g ev 
and be expressed in terms of this function v ev as 

follows: 



2E, 
A 



-(1 - 2v ev ) . 



2E K 



(1 - 2u ev ) 



(32) 
(33) 



In fact, the definition of v ev has been chosen such that 
these relations resemble the well-known expressions for 
g and k in equilibrium, where v ev has to be replaced by 
the Fermi distribution function for quasiparticles, f(E) 
(see Sec.EBJ. With the help of Eqs. ^ and the 
remaining two equations, i|28b[) and l|30|l . take the rather 
simple form 

Qod = {E ev i V 'ev } + {hod, Qod}, (34a) 
{E ev , g d} + {h odi ^ev 

}• (34b) 

Since the first of these equations is purely time-odd while 
the second one is purely time-even, we can add both equa- 
tions without any loss of information. The result can be 
written as 

j> = {£>}, (35) 
where we have introduced the new functions 



Sod 



E = E„ 



h 



od • 



(36) 



Eq. IMol) resembles very much the usual Vlasov equa- 
tion for the normal Fermi gas, which can be written as 
g = {h, g}. One just has to replace the distribution 
function g by the quasiparticle distribution function v 
and the hamiltonian h by the quasiparticle energy E. 
It should be mentioned that Eq. (|35|l or similar kinetic 
equations have already been derived in the literature sev- 
eral times. Probably for the first time it was given by 
Betbeder-Matibet and Nozieres 01 m a linearized form 
for small deviations from equilibrium. In order to be 
self-contained, we gave here our own way to arrive at Eq. 

In order to obtain a closed system of equations, Eq. 
(|35|1 must be complemented by an equation for the so-far 
unknown phase (f>. As stated above, the phase is fixed 
by the requirement that the gauge transformed gap A is 
real, i.e., ImA = 0. With the help of the relation J2jJ) 
and of the gap equation this can be rewritten as 



d 3 p 
(2nhy 



(g ev - {h ev , Qod] ~ {hod, Qev}) = . (37) 
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As we will see in a moment, this is nothing but the 
continuity equation. This observation confirms earlier 
statements in the literature that the continuity equation 
should be used for the determination of the phase [L| ■ In 
order to derive the continuity equation from Eq. I|37l) , we 
write down explicitly the poisson brackets and integrate 
by parts. In this way we obtain 

Using Eq. (|16[) and changing the integration variable ac- 
cording to p — > p + fi,V0, this can be transformed into 
the usual continuity equation, 

p(r) + V-j(r) = 0, (39) 

with 

B. Linearization around equilibrium 

From now on we will assume that the external potential 
V ex t can be written as 

Vext = Voext + Vlext , (41) 

where Vo ex t is time-independent and V\ ex t can be consid- 
ered as a small perturbation. The equilibrium quantities 
corresponding to the static potential Vo ex t will be marked 
by an index "0" , e.g., 

ivo(r,p) = /[ J Bo(r,p)] ! (42) 
to(r) = 0, (43) 

where f(E) denotes the Fermi function 

= e s/(Jr) + i ( 44 ) 

and 

E a (r , p) = ^(r, P ) + A2(r) , (45) 
P 2 

ho{r,p) = - \-V 0ex t(r) +gpo(r) - fx, (46) 



etc. Our aim is to calculate the small deviations from 
equilibrium induced by the perturbation V\ ex ti which will 
be marked by an index "1" . To that end we linearize the 
equation of motion (|35[) for the quasiparticle distribution 
function: 



i> 1 -{E ,v 1 } = f'(E ){E 1 ,E }, (47) 



where f'(Eo) = df/dEo- We also linearize the continuity 
equation i jSSjl : 



pi(r) + V ■j lv (r) - -V • Po (v)VMr) = 0, (48) 
m 

with 



j -« = /(2W^ l(r ' P) - ^ 



In order to have a closed system of equations, we must 
express E\(r, p) and yOi(r) in terms of equilibrium quan- 
tities, the perturbation Vi ex t(r), and the unknown func- 
tions vx(r, p) and </>i(r, p). Linearizing E(r, p), one ob- 
tains 



E 1 = ^-h lev + ^K l +h lod , (50) 



with 

^ie„(r,p) = V lext (r) + gpi{v) - fi<&i(r) , (51) 

h lo d(r,p) = -—p- V^(r). (52) 
m 

The most difficult part is to derive the expressions for 
pi(r) and Ai(r). We start by linearizing Eqs. 1|32|) and 



Qiev = ^-viev + - — tL^ Q [~Aq {Vlext + 9Pi - Hi) + /loA Ai] , (53) 
A 1 - 2f(E ) . 2 - 1 - 2f(E ) ~ 

Klre = -~^rV\ev 77^3 [h ^o{Vlext + 9Pl ~ Hi) + &0 A l\ + 77^ A l ' ( 54 ) 

According to Eqs. (|10fl and l|ll|) . p\ and Ai can be obtained by integrating Eqs. (|53|l and l|54|) over p. This gives a 
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coupled system of two linear equations, 



Pl (r) = p lv (r) - A{v)[V lext {v)+g Pl {v) - 7%(r)] + S(r)Ai(r) , 
Ai(r) = Ai„(r) + gB{r)[V lext (r) + g Pl (r) - 7%(r)] + [gA(r) + l]A x (r) 



(55a) 
(55b) 



r 



where the gap equation for the equilibrium case has 
been used in the derivation of the last term, and the 
following abbreviations have been introduced: 



d 3 p h (r,p) 



A(r) = A 2 (r) 



(2nh) 3 Eo(r,p 
d 3 p A (r) 



^i™(r,p) 
) 



(27rfi) 3 ^ (r,p 

d 3 p l-2f[E {r,p 



B(r) 



Ao(r) 



d 3 p 



h (r,p 



2E 3 (r,p) ' 
l-2/[£ (r,p) 
2£ 3 (r,p) 



(56) 
(57) 
(58) 
(59) 



Below we will show that the coefficient B is negligible 
compared with the coefficient A. In the limit B — > 
the two equations (|55a|l and i|55b(l are decoupled and can 
immediately be solved for p\ and Ai: 



Plv {v) ~ A{v)[V lext {v) - h^iy)} 



Ai(r) 



<?A(r) 

Ai„(r) 

gA(r) 



(60) 
(61) 



We will now calculate the coefficients A and B for the 
case that both Ao(r) and ksT are small compared with 
the local Fermi energy ,30] 



£F r 



Pf(r) 
2m 



= fJ. - V ext(r) ~5Po(r). (62) 



In this case, the relevant contributions to the integrals 
(|58|l and l|59|l come from momenta near the Fermi sur- 
face. As usual, the integrals over p can be simplified by 
transforming them into integrals over the energy vari- 
able £ = p 2 /2m — ej?(r) and approximating the den- 
sity of states by its value at the Fermi energy, i.e., 
p 2 dp « mpF(r)d^. For the coefficient A, one obtains 
in this way 



A(r) 



mp F jr) 
2ir 2 H 3 



(63) 



where the function ip describes the temperature depen- 
dence: 



cp(r) 



(64) 



One can show that ip = for T = and y> = 1 for Ao = 0. 
In all other cases, the function ip must be evaluated nu- 
merically. From its definition one can see that ip depends 
1.2 




0.4 0.6 0.8 
T/T c 

FIG. 1: Temperature dependence of the function p denned in 
Eq. fiSH . 



on r only through the dimensionless parameter T/T c (r), 
where T c (r) = 0.57A (r;T = 0)/k B is the local critical 
temperature |3l|. For illustration, the numerical result 
for ip as a function of this parameter is shown in Fig.^ 

If one applies the same method to the coefficient B, 
one obtains B = 0. This is because the integrand in Eq. 
(|S^|l is odd in £ if one neglects the energy dependence of 
the density of states. Already from this argument one 
can conclude that the coefficient B must be suppressed 
by at least one power of A /cf or T/ep. Indeed, after a 
rather lengthy and delicate analysis one finds 



B(r) 



A 



o(r) / mp F (r) 
2e F (r) V 2n 2 h 3 1 



■<P{*)] 



(65) 



This is the justification for neglecting the coefficient B 
when solving Eqs. (|55a|l and 1551 



Finally, let us put everything together and give a con- 
cise summary of the system of equations which has to 
be solved. First of all, there is the equation of motion 
(|47|l for the quasiparticle distribution function. After 
the Poisson bracket on the r.h.s. has been written down 
explicitly, it can be transformed into 
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!>1 - {E Q ,Vi} 



f'(E ) 



.Vlext + gPlv - fr<t>l , Ao „ Ao{Viext + 9Plv ~ ft(/>l) , ^0 
J 



pV — + ^ P ' V TT^i + ^ p ' v 



A Ai„ 
5^4 



--^(P- V) 2 <j> 1 -h(^V(V 0e:c t+gpo) + ^VAo) ■ Vfa . (66) 



The second equation is the continuity equation 

Pl „(r) - A{r)[V lext (r) - hMr)] 
l+gA(r) 



V-Ur)--V-p„(rWi(r) = 0. 

m 



(67) 



The definitions of pi„, Ai^, and j lv in terms of v\ are 
given by Eqs. (JSSJl, l(57)l. and (gSJ. 



C. Limiting cases 

We are now going to check that our equations repro- 
duce superfiuid hydrodynamics and the Vlasov equation 
in the cases T = and T > T c , respectively. In the limit 
of zero temperature, Eq. I|66|l becomes extremely simple 
since f(E) = and therefore the r.h.s. of Eq. (|66|) van- 
ishes identically. The corresponding solution is of course 
v\ = |32|, which implies p\ v = A^ = j ly = 0. As a 
consequence, the continuity equation (|67f) reduces to 



Viextjr) - h'(f>i(v) h_ 
2ir 2 H 3 m 



V-po(r)V0i(r)=O. (68) 



mp F (r) 



Here we have used the explicit expression for A(r) and 
the fact that ip = at zero temperature. 

How does Eq. (|68[) compare to superfiuid hydrodynam- 
ics? The continuity and Euler equations of superfiuid 
hydrodynamics can be written as |5j]: 



v(r) = -V 



p(r) + V • p(r)v(r) = , 
v 2 (r) Vext(r) , Mioc(r) 



m 



m 



(69) 
(70) 



where v(r) denotes the velocity field and pioc( r ) is the 
local chemical potential, which in the BCS phase (A <C 
€f) is related to the density p(r) by the Thomas-Fermi 
relation 



p 2 F ( r ) 

2m 



9p(r) 



(71) 



with 



Mr) = ft[67r 2 p(r)] 1/3 . (72) 
Writing the irrotational velocity field in the form 



v(r) = V0(r) 

to 



(73) 



and linearizing Eqs. I|69|) and (|70|) around equilibrium, 
one obtains 



Pi(r)--V-po(r)V&(r)=0, 
m 



fi0(r) - 7i« 



t (r) + ( 



2ir 2 h 3 



(74) 
(75) 



>mj)F(r) 

Solving Eq. I|75[) for p\ and inserting the result into Eq. 
I|74jl. one reproduces exactly Eq. (|68|l . This can be seen 
as an alternative to the recent derivation of superfiuid 
hydrodynamics from the underlying microscopic theory 
in Ref. 

The analysis of the other limit, T > T c , is more diffi- 
cult. In this limit, the gap Ao vanishes and consequently 

£ (r,p) = Mr )P )|, (76) 
viev(r,P) = sgn[p-p F (r)]£i e „(r,p) . (77) 

In addition, one has <^(r) = 1, A(r) = 0, pi(r) = pi„(r), 
and Ai(r) = Ai„(r) = 0. Using these relations, and 
considering separately the two cases p < pf (i.e., ho < 0) 
and p > pf (i.e., ho > 0), one can convince oneself that 
Eqs. iJHSJ and (JSTJ) reduce to 



Qi - {ho, Qi} 

h 



f'(ho 



m 



-(P-V) 



( - p ■ V{Vi ext + gpi - U4>i) 
h[V(Vo ext +gpo)} ■ V^i 



(78) 



and 



Pi(r) + V .j lv (r) - -V • p (r)V^i(r) = 0. (79) 

TO 

As we will see in a moment, these two equations are not 
independent of each other. Hence, they do not allow 
to determine £>i(r, p) and (/>i(r) in a unique way. This 
is in fact very reasonable since the condition Im A = 
fixing the phase <fr becomes meaningless above T c , where 
A = 0, and therefore the function <fi should be completely 
arbitrary in this case. The relevant physical quantity, 
which of course should be unique, is f?i(r, p). Linearizing 
Eq. (|Po)l and using £>o( r ,p) = f[ho(r,p)], we can express 
g>i(r, p) in terms of Qi(r, p) as follows: 



0i(r,p) = £>i(r,p) /'[/i (r,p)]p- V<; 



(80) 
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If we insert this into Eq. JTSJl, all terms containing the 
phase 0i drop out, and we are left with 

Qi - {h , qi} = /'(M- ' V(VW + g Pl ) . (81) 
m 

This is nothing but the linearized form of the Vlasov 
equation, 

Qi - {ho, Qi} = {hi, qq} , (82) 

with hi = Vi ex t + gp±. It remains to check that the con- 
tinuity equation (|79|l is satisfied for arbitrary functions 
4>i, if Qi fulfills Eq. lj%T)> . To that end, we multiply Eq. 
(f£T|) by p and integrate over p, which leads to the usual 
continuity equation 

p 1 (r) + V-j 1 (r) = 0. (83) 

With the help of Eq. (|8C)|> the current j x can be written 
as 

h(r)=hu(r)- -po(r)V<Mr). (84) 
m 

Inserting this into Eq. (|83|l . we indeed recover Eq. (|79|l . 
Since we did not make any assumptions about the func- 
tion </>! (r) , we conclude that it is completely arbitrary, as 
it should be. 



III. SIMPLE EXAMPLES 

A. Sound wave in a uniform system 

In this subsection we are considering a particularly 
simple excitation, namely a sound wave traveling through 



a uniform medium. This case has already been studied 
by Leggett 0] many years ago (except for the numerical 
evaluation of the integrals) by using the standard tech- 
niques of normal and anomalous Green's functions. The 
purpose of the present subsection is therefore to check 
that our apparently very complicated equations H66(l and 
()67J) correctly interpolate between the limits of zero and 
critical temperature. 

Since the medium is assumed to be uniform, the equi- 
librium quantities do not depend on r. We consider an 
excitation operator of the form 

V lext (r;t) = V lext e ik - r - i ^. (85) 

As usual, in order to ensure that the perturbation van- 
ishes for t — > — oo, one can assume that oj has an in- 
finitesimal positive imaginary part. From translational 
invariance it is clear that all quantities describing the 
deviations from equilibrium will also have the form of a 
plane wave, with the same wave vector k and frequency 
u> as the excitation. Like V\ ex t, the amplitudes will be 
marked by a hat over the corresponding symbol. Con- 
cerning the phase </>i, it turns out to be convenient to 
parametrize it in the form 

<f> 1 {r;t) = 4> l -J*- r - iut . (86) 

UJ 

The Poisson bracket on the l.h.s. of Eq. (ffjfjfl now becomes 

rp i • ho P - k ik-r-iwt / Q n\ 

{E , v\) = -i— v x e , (87) 

Ep Tfl 

and Eq. (|66(l can easily be solved for v\ : 



_ ,, - P • k h / hp V lext - tkjyx + gp lv Ao Ai^ p ■ k ? \ 
^ °' muj E \E 1 + gA Ep gA muj J 

hpTp-'k 

Ep muj 



Vi — . (o»J 

hp p • k 



Of course, the quantities p\ v and K\ v on the r.h.s. de- 
pend themselves on v\. Therefore the next step consists 
in inserting this expression for £>i into Eqs. H5(jfl and (|57|l . 
The integrals over the angle between p and k can be eval- 
uated in closed form. For the remaining integrals over p, 
we will again exploit the fact that the gap and the tem- 
perature are much smaller than the Fermi energy, as we 
did already in Sec. Ill Bl We thus replace p 2 dp by mpFd£, 
and in the integrand we replace p by pf, except for hp 
and Ep, which must be replaced by £ a nd ^/ g 2 + Aq, 
respectively. Like the coefficient B in Sec. Ill Bl the inte- 



grals which lead to the coupling between the equations 
for p\ v and Ai„ are zero within this approximation, i.e., 
they are of higher order in A/ep or T/ep and can be 
neglected. The resulting equation for piv reads 



Plv 



mpp / V lext - hipi + gp lv 



2n 2 h 3 



1 + gA 



h{s) + h<j) 1 I (s) 



(89) 
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Here we have introduced the abbreviation 



In(s) 



E 



sE 

T 



arctanh 



r (90) 

where E = ^/ £ 2 + A 2 ,, and s denotes the dimensionless 
ratio of the sound velocity c = w/k and the Fermi velocity 
vf = Pp/m, 



c muj 
vp ppk 



(91) 



Although not marked explicitly, I n (s) depends not only 
on s but also on the ratio T/T c (analogously to the func- 
tion <p). Note that the integrals I n (s) have a branch cut 
along the real axis from s = — 1 to s = 1. The infinites- 
imal imaginary part of u>, i.e., of s, fixes the sign of the 
imaginary part of I n (s). 

Until now we have one equation for two unknown quan- 
tities, p\ v and 4>i. The second equation can be obtained 
from the continuity equation (|67|) . It is evident that the 
current flows in longitudinal direction, such that it 
can be written in the form 



(92) 



We will now express ji u in terms of Viext ^ 0ij and p\ v 
by inserting Eq. (|88|l into Eq. I|49|) . The integration over 
p is done as explained above for the case of /5i„, and the 
result reads 



Jlu = 



mCpF ( Viext - fr<t> + gPlv 



lir 2 ^ \ 



1 + gA 



h<> J A s))-^. (93) 



mc 



In the last term, we have introduced the "normal density" 
of the system, p n , which is given by 



Po - Ps 



-po J dZf'(E ) 



(94) 



Correspondingly, p s is the "superfluid density" . Note 
that the ratios p n /po and p s /po depend only on one pa- 
rameter, namely T/T c . The numerical results for p n /po 
and p s /po as functions of T/T c are shown in Fig. [2 

Inserting Eq. jHSjl into the continuity equation ljfl7jl . 
one obtains the second equation which is needed for de- 
termining p\ v and 4>x '■ 



plu - A(Vl e xt - fl4>l ) . rriPF ( Vi ext - fr<fi + Qplv T , , fe t T , A Paf>4> 



( Viext ~ M> + gp lv I \ PstKp! 

I — — — J (a) + ^i7_ 2 (s)) r=°- 95 ) 

V 1 + a A / mc z 



1 + gA 2tt 2 ?i 3 V 1 + gA 

In principle we could now solve Eqs. (|89|l and (|95|) for the two unknown variables p\ v and 4>i. However, it is more 
transparent to use the amplitude of the total density oscillations, pi, as variable instead of the auxiliary quantity p\ v . 
Expressing p\ v in terms of pi with the help of Eq. I|60[) . we rewrite Eqs. I|89|) and (|95|l as 



g|§[l -<p + h{ S )])pi - ^[1 -<p + I 2 (s) - I (s)]4l = [1 -V + h{s)]Vi ext , (96a) 

3s 2 po- 



, + 2^F /0(S) > X + WV-^ /o(s) - 3^ J^ 1 = ~2^¥ Io{s)Vlext ■ (96b) 



r 



It is straight-forward to solve this 2x2 system of equa- 
tions for pi. Let us introduce the response function n, 
defined such that 



Pi 



Tripp 
2ir 2 h 3 



U(s;T/T c ;k F a)Vi. 



(97) 



The first term has been factored out in order to make 
n dimensionless. From the system of equations <|96[) one 
can see that IT is a function of s and the two parameters 
T/T c and 



Lo = 



gmp F 



(98) 



The explicit expression for II can most conveniently been 



expressed in the form 



II(s; T/T c ; k F a) = 



n (s;T/T c ) 



2kFd 

1 U (s;T/T c 



(99) 



where IIo is the response function in the limit kpa — > 0: 



n = 



3s 2 po 



}_Ps 

3s 2 po 



(100) 



1-77^7 <P + h + 1-2 ~ 2J 



Note that these expressions coincide exactly with the 
quantum mechanical result in the long-wavelength and 
low-frequency limit as given by Eqs. (68) and (69) of Ref. 



10 




FIG. 2: Temperature dependence of p n /po (solid line) and 
p s /po (dashed line). 
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T/T c = 0.99 ----- 
\ T/T c =1 
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FIG. 3: Excitation spectrum — Imll of a uniform medium as 
function of the reduced sound velocity s — c/vf for different 
temperatures. 



12] . In order to see this, it is sufficient to observe that 
after integration over the solid angle the quantities a, £, 
and r\ defined in Eq. (65) of Ref. ^3 can be expressed 
in terms of our integrals as a = (1 — <p + 72 — 7o)/2, 



C - [p s /(3p ) + s 2 (I Q - J_ 2 )]/2, and V 



-Iq. In our 



case of a pure s-wave interaction, the Landau parame- 
ters in Eq. (68) of Ref. [EJ are given by Fq = 2kpa/Tr 
and F\ = 0. Then the quantities K\ and Q of Ref. [IJ 
correspond to our IT and n , respectively. As stated in 
Ref. [lj , the long- wavelength and low- frequency limit is 
valid if hio^vphk <C A. Since our semiclassical result co- 
incides with this limit of the quantum mechanical result, 
we conclude that this is the condition for the validity of 
our semiclassical theory. A calculation of the response 
function beyond the long-wavelength and low-frequency 
limit can be found in Ref. |23|. 

The excitation spectrum of the system is characterized 
by the imaginary part of II, which is plotted in Fig. [31 for 
kpa = —0.25 and several temperatures between 0.8 T c 
and T c . One can see that at 0.8 T c the excitation spec- 
trum exhibits a peak near s ~ 0.5 which becomes broader 
and finally disappears when the temperature approaches 
T c . 

In order to interprete this behavior, let us again con- 
sider the two limits T — > and T — > T c . In the zero- 
temperature case, all integrals containing the term /'(Eq) 



in the integrand vanish, i.e., <p — p n /po — 7 n (s) = 0, and 
the response function reduces to 

U(s;0;k F a)= 1 . (101) 

3s z — 1 — 2kpa/TT 

This means that the excitation spectrum is a 6 function 
at 



1 2kpa 

- A — 

3 



3tt 



(102) 



corresponding to the hydrodynamic speed of sound. 

In the other limit, T — > T c , one has if = p n /pa = E 
The integrals I n (s) reduce to 

IJs; T/T c > 1) = 1 - s arctanh - , (103) 

s 

independent of n, since the factors £/E in the integrand 
of Eq. H90f) can be replaced by 1. As a consequence, the 
two equations of the system (|96|) become identical and 

the coefficients in front of <f>\ vanish, in accordance with 
the more general arguments of Sec. Ill CI Solving for pi 
gives 



IL{s;T/T c > I- k F a) 



1 



— ^1 — s arctanh — 
2kpa 



zpa / 

7T \ 



1 



s arctanh — 
s 



(104) 

in agreement with the usual result of Landau's Fermi- 
liquid theory for the case of a pure s-wave interaction. 
If the interaction was repulsive (a > 0), Eq. (|1()4|) would 
have a pole at s > 1, corresponding to the propagation 
of zero sound. However, here we are considering the case 
of an attractive interaction, where zero sound does not 
exist. Instead there is a continuous spectrum of particle- 
hole excitations ranging from s = to s = 1. 

Our numerical results shown in Fig. [3] can be inter- 
preted as follows. At zero temperature, there exists a 
collective hydrodynamic sound, which is undamped (at 
least within the present theoretical treatment). As the 
temperature increases, a normal component consisting 
of thermally excited quasiparticles builds up. However, 
at temperatures where p n is already considerably differ- 
ent from zero, the hydrodynamic sound is still practi- 
cally undamped. The reason for this is that all thermally 
excited quasiparticles contribute equally to p n , whereas 
only those quasiparticles whose velocity dE/dp « vp^/E 
is at least equal to the sound velocity c contribute to 
the Landau damping. At sufficiently high temperature, 
the Landau damping becomes very strong and the hy- 
drodynamic sound ceases to exist. What remains is a 
continuum of particle-hole excitations, and the interac- 
tion manifests itself only in the rounded edge near s = 1 . 



B. Quadrupole mode in a spherical trap 

Our main motivation for developing the present semi- 
classical approach was to apply it to the case of trapped 
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atomic Fermi gases. The simplest example which comes 
to our mind is the quadrupole oscillation of a Fermi gas 
in a spherical trap. Even in this case, the r dependence of 
the equilibrium quantities induced by the trap potential 
makes our equations very complicated. Since in this first 
investigation we are interested in problems which can be 
solved analytically, we will apply two additional simpli- 
fying approximations, which allow us to obtain explicit 
solutions. A numerical method for solving our equations 
without additional approximations will be proposed at 
the end of this subsection. 

Let us start with the linearized equation (|66|) . which 
has the form 

!> 1 (r,p;t)-{£; (r,p),i/ 1 (r,p;t)} = F(r ) p;t). (105) 

For its solution we adopt the Green function method used 
in Ref. [23| to solve the linearized Vlasov equation for 
nuclear giant resonances, which is formally very similar to 
our problem. The starting point is to write the solution 
of Eq. (|105|) in the form 



H(r,p;t) = J dt' f d 3 r' [ d 3 p'G(r,p,r',p';t-t') 

xF(r',p';t'), (106) 

where G is the Green function of the differential operator 
on the l.h.s. of Eq. (|105fl . satisfying 



A_ v ( 9E ° 9 9E ° 9 

dt ^ V On dpi dpi dn 



G(r,p,r',p';t) 



8(t)5(r-r')5(p-p'). (107) 



Denoting by R(r, p; t) and P(r, p; t) the solutions of the 
classical equations of motion 



Ri 



0E Q (R,P) 

dPj 



Pi = (108) 



satisfying the initial conditions R(r, p; 0) = r and 
P(r, p; 0) = p, one can show that 

G(r, p, r', p'; t) = 0(t)S[r - R(r', p'; t)]8[p - P(r', p'; t)] . 

(109) 

fulfills Eq. I|107|) . Due to time-reversal symmetry and 
Liouvillc's theorem, this Green function can be rewritten 



G(r, p, r', p'; t) = 0(i)<J[r'-R(r, -p; i)]d[p'+P(r, -p; t)} . 

(HO) 

The latter form renders the phase-space integrals in Eq. 
(|1U6[) trivial. Changing the time integration variable ac- 
cording to r = t — t', one obtains 

/>OC 

M(r,p;t)= / drF[R(r,-p;r),-P(r,-p;r);i-r]. 
Jo 

(HI) 

In the case of a harmonic perturbation, 



V lext (r;t) = V lext {v)e~ 



the time dependence of v\ as well as the explicit time 
dependence of F will be harmonic, too. We will denote 
the amplitudes by a hat over the corresponding symbols. 
Multiplying Eq. I|lll|) by e luJt , one finds that v\ is given 
by the Fourier integral 



*>i(r,p) 



dr e toT J'[a(r, -p; r), -P(r, -p; r) 



(113) 

For the purpose of illustration we want to discuss a 
simple case in which the classical trajectories are ana- 
lytically known. We make two approximations: First, 
we replace the r-dependent gap A (r) by a constant A . 
This approximation implies that v\ ev is odd in £ and 
Ai„ can be neglected, as it was the case in the preced- 
ing subsection. Second, we will neglect effects from the 
Hartree mean-field as well in the equilibrium state as in 
the deviations from equilibrium. The second approxima- 
tion, which is by far not as unrealistic as the first one, 
amounts to neglecting all gpo and gp\ v terms and replac- 
ing the denominators 1 + g A in Eq. (|<j fc>|) by 1. The trap 
potential is assumed to be a spherical harmonic oscilla- 
tor, 



1 



Voext = mfl r 



2^2 



(114) 



It is evident that the equations of motion (|108|l conserve 
Eq. However, if Ao is a constant, this implies that ho is 
conserved, too, and the solutions of Eqs. I|108fl are closely 
related to the those of the ordinary harmonic oscillator. 
Indeed, it is straight-forward to show that the trajectories 
are given by 



R(r,p;f) = rcos 



ho fit 

Eq ' mfl 



p ho^lt 
sin ■ 



En 



h Q,t . h Qt 

P(r, p; t) = p cos m\LY sin — — — 

Eq Eq 



(115a) 
(115b) 



Since ho and Eo are constants of the motion, they can 
likewise be evaluated at (r,p) or (R, P). 

Due to our approximations, the function F [given by 
the r.h.s. of Eq. (|66|l ] reduces to 



F 



-f(Eo) 



hi p 



El m 



V(Vlesi + ihLJ(f)i 

-.2 ho 



(112) 



v tyi-sr^T-vtufn . (lie) 

As excitation we choose the quadrupole operator 

t> lelt (r)=amO 2 (r®r) 20 , (117) 
where, explicitly, 

I \ 1 IOn\ 2V Z W Z ~ V X W X - VyWy 

(v(giw) 20 = > (l^lv\20)v fj,w u = — . 

— ' v6 

(118) 

The prefactor in Eq. (|117|l has been chosen such that 
the coefficient a is dimensionless. Due to the spherical 
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symmetry of the trap, the angular dependence of (f>i must 
be of the same quadrupolar form as that of V\ ex t, but the 
radial dependence could in principle be different. Here 
we make the ansatz that <f>\ is proportional to Vi ext , i.e., 

<£i(r) = /?^(r®r) 2 o, (119) 

and we will show afterwards that with this ansatz for <p 
the continuity equation can be satisfied by an appropriate 



choice of the coefficient (3. Quadratic ansatze like Eq. 
(|119fl . corresponding to a superfluid velocity field which 
is linear in the coordinates, have frequently been used 
(see, e.g., Ref. 0) for the calculation of the frequencies of 
collective modes in the limit of superfluid hydrodynamics 
(T = 0). 

Inserting Eqs. 1117(1 and 1(119(1 into Eq. 1(116(1 and using 
the explicit form of the trajectories, Eq. ((115(1 . we can 
evaluate the Fourier integral in Eq. 1(1131) , with the result 



v\ = 



= -2/'(£ ) 



(P ® P)20 



hi 



-mfr(r®r) 20 — 

£j 



ho fl Eq Eq 



El 



r 



hi 



n(r®p) 20 -^ 



Aq 

Eq 




IV 2 

ft 2 " 


hi 

4— 

El J 



(120) 



Now we have to calculate the corresponding current j and 
density oscillations p\ u . As detailed in the preceding sub- 
section, this is accomplished by integrating pvi/m and 
z)i, respectively, over p. Replacing in the integrals p 2 dp 
by mp F (r)d£, p by pj?(r), h by f, and E by a/£ 2 
we obtain 



A 2 



hu( T ) = Po - 4J 22 (z)] - iazl 20 (z)) V(r ® r) 20 , 

(121) 

i,r[l ~ l " <r) [aho(z) ~ ipzl22(z)](r ® r) 20 , 

(122) 



with the abbreviations 



uJ 

n 



J« 0) 



fA 1 

z 2_ 4 £2/£2 ' 



(123) 
(124) 



Aq. From its definition it is evident 



where E = 

that /40 = I20 — ^22, such that it is sufficient to evaluate 
two of these integrals numerically. The functions /y (z) 
have a branch cut along the real axis from z = —2 to 
z = 2. Remember that ui, and therefore also z, is assumed 
to have an infinitesimal positive imaginary part, fixing 
the sign of the imaginary part of hj(z). 

As stated above, the coefficient (3 must be determined 
by the continuity equation 1(67(1 . Due to our approx- 
imation to neglect the Hartree field, the denominator 
1 + gA(r) in the first term of Eq. 1(67(1 can be replaced by 
1, and the Fermi momentum can be given in closed form: 



Pf(t) = 



/! — —mVl 2 r 2 



(125) 



Inserting the results for ] lv and p\ v into the continuity 
equation, one finds that the ansatz l(119|) indeed allows 



to satisfy the continuity equation, and the corresponding 
solution for the coefficient (3 reads 



/3 — iza 



1-^ + 2/22(2) 



(l-<p)(z 2 -2) + 2I 22 (z)(z 2 -4) 



(126) 



This expression can be used to obtain p\ v . Here we will 
immediately give the result for the amplitude of the total 
density oscillations, i.e., p\ = p\ v — A(Vi ext + itkofii), 
which we write in the form 

Mr)=a m2 "y\ r®r) 20 n(z) (127) 



with 



U(z) = ho(z) + 



[l-<p + 2h 2 (z)][l-<p + 4/ 22 (2)] 



(l-^)(z 2 _2) + 2/ 22 (z)(z 2 -4) • 

(128) 

Before discussing numerical results, let us again study 
the two extreme cases T = and T > T c . In the zero- 
temperature limit, all integrals tp and /y are zero, and 
hence the response function becomes 



n(z;T/T c -0) 



1 



(129) 



z 2 - 2 ' 

i.e., it has a single pole at the hydrodynamical frequency 



V2~n. 



(130) 



In the case T > T c , i.e., in the normal phase, we have 
(p = 1, and in the definition 1124(1 we can replace Ao and 
Eq by and £, respectively, such that we obtain 



1 



/20 (2; T/T c > 1) = 



h2{z-T/T c > 1) =0. 



(131a) 
(131b) 
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FIG. 4: Quadrupole excitation spectrum — Imll of a harmon- 
ically trapped gas as a function of the excitation frequency (in 
units of the trap frequency) for different temperatures. The 
spatial dependence of the gap as well as the Hartree mean- 
field have been neglected. 



Thus the response function reduces to 

1 



U(z;T/T c >l) 



z 2 -4 



(132) 



Like in the zero-temperature case, we have a single pole, 
but now at a frequency which is higher by a factor of y/2. 
The reason for the difference of the two frequencies is as 
follows. In the superfluid phase, the momentum distri- 
bution stays spherical during the oscillation. In contrast 
to this, in the normal phase, the momentum distribution 
is deformed in the opposite direction as the density in 
coordinate space. This deformation of the Fermi sphere 
costs kinetic energy, which increases the restoring force 
and thereby the frequency of the oscillation. 

At intermediate temperatures < T < T c , the exci- 
tation spectrum is continuous and it is characterized by 
the imaginary part of II (z), which is shown in Fig. 0] for a 
set of temperatures between 0.5 T c and T c . At 0.5 T c , the 
spectrum exhibits a sharp peak at the hydrodynamic fre- 
quency z = V2, the weak broadening being due to Lan- 
dau damping. With increasing temperature, the Landau 
damping becomes more important, and at the same time 
the centroid of the distribution moves to higher frequen- 
cies. Above 0.8 T c , however, the width of the peak does 
not increase any more with temperature, but it decreases. 
Finally, when the temperature approaches T c , the peak 
becomes again very sharp and, not surprisingly, it lies at 
the frequency z — 2 predicted by the Vlasov equation. 

One might ask the question whether the approxima- 
tions made in this subsection are justified or not. Let 
us therefore compare our results with those of a QRPA 
(quasiparticle random-phase approximation) calculation 
pj, where, apart from the mean-field approximation lead- 
ing to the Bogoliubov-de Gennes equations, no approxi- 
mations are made. Qualitatively our scmiclassical results 
show the main features of this quantum-mechanical cal- 
culation: the hydrodynamic mode at zero temperature, 
its damping at intermediate temperatures, and the sub- 
sequent reappearance of an undamped collective mode 
with a higher frequency in the normal phase. That our 



frequency in the normal phase is exactly equal to z = 2 is 
a consequence of neglecting the Hartree mean field. How- 
ever, in the range of validity of our theory (&fM <C 1), 
the Hartree mean field cannot shift the frequency very 
much (in Ref. Q, e.g., the frequency is shifted from 2 to 
« 2.2). We therefore believe that this effect is not very 
important. More problematic is the constant-gap approx- 
imation which we needed for the analytical solution of the 
equations of motion (|108fl . Because of this approxima- 
tion, there are no quasiparticlcs having energies below 
A , and as a consequence, the Landau damping sets in 
at rather high temperatures. In the full calculation, how- 
ever, the lowest-lying quasiparticles are those whose wave 
functions are localized near the surface, where the gap is 
small, and which have much smaller energies than the 
central value of the gap. Therefore within the full cal- 
culation the Landau damping is already quite important 
at very low temperatures. In the semiclassical formal- 
ism these low-lying quasiparticles can be understood as 
quasiparticles bouncing back and forth between the po- 
tential wells created by the trap potential and the spa- 
tially varying gap ( Andreev reflection) [23, [25| . The in- 
clusion of this effect would require a numerical solution 
of the equations of motion l|108|) . 

This leads us to a possible numerical method for solv- 
ing even the original (i.e., not linearized) kinetic equation 
(|35|l . In nuclear physics, the Vlasov equation (usually 
complemented by a collision term) is routinely solved by 
the so-called test-particle method, e.g., in order to sim- 
ulate heavy- ion collisions |2g. Recently this method has 
also been applied to the solution of the Vlasov equation 
with collision term for trapped atomic Fermi gases |27| 
and of a Vlasov-like equation for trapped fermion-boson 
mixtures |2^| , and it seems to be straight-forward to gen- 
eralize it to our case. The basic idea of the method is as 
follows. Instead of calculating the time evolution of the 
continuous quasiparticle distribution function v(r, p;t), 
one can use a finite number of "test-quasiparticles" and 
follow their motion in phase-space by solving numerically 
the equations of motion 



Ri — 



d£(R,P) 

dPi 



Pi = - 



dE(R, P) 



(133) 



for all test-quasiparticles simultaneously. Of course, the 
quasiparticle energy E(r,p;t) contains the mean-fields 
gp(r;t) and A(r;i), which must be calculated at each 
time step from the actual quasiparticle distribution. At 
the same time, the phase </> must be calculated at each 
time step from the continuity equation. This seems to be 
a tractable task, which will be adressed in a subsequent 
publication. 



IV. SUMMARY AND CONCLUSIONS 

In the first part of the present article, we derived a 
set of semiclassical equations describing the dynamics of 
a collisionless superfluid Fermi gas by taking the h — » 
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limit of the TDHFB or time-dependent Bogoliubov-de 
Gennes equations. In the limits of zero and critical tem- 
perature, these equations reproduce superfiuid hydrody- 
namics and the Vlasov equation, respectively. At inter- 
mediate temperatures, there is a complicated interplay 
between the dynamics of the superfiuid component of the 
system, governed by the function <p(r) which is related to 
the phase of the gap A(r), and the dynamics of the nor- 
mal component, which is described by the quasiparticle 
distribution function v(r, p). The dynamical equation 
for v formally corresponds to the usual Vlasov equation 
with g and h replaced by v and E, respectively, while 
the function <fi is determined by the continuity equation. 
The latter point can be seen most easily in the linearized 
version of the equations for small deviations from equi- 
librium. 

In the second part, we gave an illustration of our equa- 
tions by applying them to two simple cases, where ana- 
lytical solutions could be found. The first example we 
studied was a sound wave traveling in a uniform system. 
In this case we could reproduce the usual hydrodynamic 
speed of sound at zero temperature. At non-zero tem- 
peratures below T c , the sound wave suffers strong Lan- 
dau damping because of its coupling to thermally excited 
quasiparticles. For T — > T c the excitation spectrum con- 
tinuously goes over into that of the usual particle-hole 
continuum (with RPA corrections) which is found above 
T c . 

The second example was the quadrupole mode of a 
Fermi gas in a spherical trap. Applying the approxima- 
tion of a constant gap and neglecting the Hartree field, 
we were able to solve the linearized quasiparticle kinetic 
equation exactly also for this case. We could qualitatively 
reproduce the most important results of quantum me- 
chanical QRPA calculations: At zero temperature, there 
is an undamped collective mode at the hydrodynamic fre- 
quency lu — V2fi, which becomes strongly damped at low 
temperatures (0 < T <C T c ). At a certain temperature 
the damping rate reaches a maximum, and above that 
temperature it decreases until at T = T c an undamped 
collective mode reappears at the frequency predicted by 
the Vlasov equation, which is higher than the hydrody- 
namic frequency. However, quantitatively the agreement 
with the QRPA calculation is not yet satisfactory, essen- 
tially because we replaced the gap by an r-independent 
constant. We suggested to use the test- particle method 
for the numerical solution of the equations in the case of a 
spatially varying gap A(r), which at the same time would 



allow to include the Hartree field and to treat strong de- 
viations from equilibrium, like the expansion of the gas 
after the trap is switched off. 

As long as the gas is close to equilibrium, i.e., as long 
as the Fermi surface can be regarded as spherical, the 
effect of collisions is very small due to Pauli blocking of 
the final states. However, in the case of strong devia- 
tions from equilibrium, like during the expansion of the 
cloud when the trap is switched off, the deformation of 
the Fermi sphere can lead to rather important collisional 
effects 29]. Therefore, in this case it is necessary to in- 
clude the collision term into the theory. This is an in- 
teresting problem which should be addressed in a future 
investigation. For normal-fluid trapped Fermi gases there 
exist already some calculations which take the collision 
term into account [53, . The more complicated case of 
paired Fermi systems with collisions has been considered, 
e.g., in the context of superfiuid 3 He [l4llT5|. 

Since our equations were obtained as the K — *■ limit of 
the TDHFB equations, only the leading gradient terms 
are included. This can be seen, e.g., in our results for 
the sound wave, which in fact are the long-wavelength 
and low-frequency limit of the quantum mechanical re- 
sult which can be obtained by diagrammatic techniques. 
For the system in a trap this means in particular that 
the gradients of the trapping potential, which are propor- 
tional to the trap frequency f2, must not be too strong. 
But what is "too strong"? In the normal phase (A = 0) 
the h — > limit works extremely well if Ml <§C [i, which 
is always the case in the experimental situations. In the 
superfiuid phase (A ^ 0) the relevant condition reads 
HQ <C A [3, j which is much more difficult to satisfy, 
especially for the radial trap frequency which is usually 
much larger than the axial one. In spite of this limitation, 
we believe that the semiclassical approach presented here 
will be useful, since at present no quantum mechanical 
calculation is able to describe the dynamics of systems 
with more than « 10 atoms, especially in the case of 
deformed traps. 
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